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Geometric  Data  Analysis: 

An  Interactive  Graphics  Program 
for  Shape  Comparison 

by 

Andrew  F.  Siegel 
ABSTRACT 

Two  shapes,  each  consisting  of  n  homologous  points, 
can  be  rotated,  scaled,  and  translated  to  obtain  a  close  fit 
to  each  other  by  several  methods.  An  interactive  graphical 
computer  program  is  presented  here  that  implements  two  methods 
least  squares  and  repeated  medians,  a  robust  method.  Examples 
are  given  and  the  use  of  the  system  is  discussed. 


1.  INTRODUCTION 


The  quantitative  study  of  shape  and  form  has  been  considered 
by  many  authors  since  the  fundamental  descriptive  work  of 
Thompson  ( 1 0 1 7 ) .  The  method  of  least  squares  was  used  by  Sneath 
(1967)  to  solve  the  problem  of  finding  a  common  location  and 
orientation  of  two  shapes  in  order  to  compare  their  similarities 
and  differences.  A  general  discussion  of  the  numerical  aspects 
of  such  least  squares  calculations  may  be  found  in  Huber  (1980). 
There  are  situations  in  which  robust  methods  such  as  the  repeated 
median  technique  (Siegel,  1980)  are  far  superior  to  least 
squares,  especially  in  the  detection  of  localized  shape  differ¬ 
ences  (Siegel  and  Benson,  1979). 

Least  squares  and  robust  methods  are  both  considered  here, 
and  an  interactive  computer  program  written  in  FORTRAN  with 
graphics  capabilities  is  presented.  Section  2  outlines  the 
mathematical  methods  involved  in  each  fitting  process;  further 
details  can  be  found  in  the  references.  Two  examples  are  provided 
in  Section  3  to  illustrate  the  use  of  the  program  and  to  show 
the  kinds  of  graphic  results  that  can  be  obtained.  Section  4 
outlines  the  work  involved  in  setting  up  the  program  to  run  on 
a  computer  installation  and  provides  an  overall  description  of 
how  to  use  the  system.  The  interactive  commands  are  described 
in  detail  in  Section  5,  and  the  program  is  listed  in  the  Appendix. 
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2.  METHODS 


* 

I 


Each  of  the  two  shapes  is  represented  by  a  sequence  of  n 
points  in  two  dimensions.  Homology  is  assumed,  so  that  the  ith 
point  of  shape  1  corresponds  to  the  ith  point  of  shape  2  because, 
for  example,  this  might  be  the  location  of  the  base  of  the 
skull  in  each  specimen.  Establishment  of  homology  can  be  straight¬ 
forward  in  some  situations,  but  is  difficult  in  others.  We  will 
assume  it  has  been  done,  even  if  only  tentatively,  and  will  also 
assume  that  any  needed  reflection  has  also  already  been  done.  Thus 


our  data  are 

Shape  1 

Shape  2 

(Xj.yj) 

(u1,v1) 

(x2,y2) 

(u2,v2) 

<V*n> 

(»«•»„) 

Holding  shape  1  fixed,  we  transform  shape  2  by  a  rotation 
angle  0  ,  a  magnification  factor  p  ,  and  a  translation  (a,b)  . 
The  transformed  coordinates  of  shape  2  are  therefore 


cos(0)  -sin(e) 
+  \sin(e)  cos(S) 


(2.1) 
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(2.7) 
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The  repeated  median  procedure  estimates  the  magnification 


(p) 

and 

the  rotation 

(9)  separately. 

using  a  double  median 

for 

each . 

These  estimates  are 

-  _  Median 
p  ~  l<i<n  < 

f Medi an 
(  l<j<n 

V,  j*i 

(2.8) 

V 

s  _  Median 

0  "  1<i<n  < 

/^Medi  an 
[  1  <  j  <n 

l 

-5 

(2.9) 

where  p .  ■  and  0-.  are  the  magnification  and  rotation  para- 

I  J  '  J 

,Sii  meters  computed  using  only  points  i  and  j  of  each  shape. 

The  transformation  factors  can  be  estimated  using  a  single  median 

i  ' 

!  a  =  {xi-5Cu1cos(e)-v1s1n(5)]>  (2.10) 

b  =  {yrp[uisin(0)+v1cos(e)]}  (2.11) 

**  The  match  obtained  using  repeated  medians,  unlike  least 

squares,  will  not  be  led  astray  by  a  few  atypical  points  and 
localized  changes  will  be  more  easily  identified.  In  general, 

r» 

if  more  than  (n+l)/2  (i.e.  just  over  half)  of  the  points  can  be 
made  to  match  closely,  then  the  repeated  median  method  will 
produce  a  transformation  that  does  match  them  closely. 


I 
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3.  EXAMPLES 

Two  examples  are  given  here  in  order  to  illustrate  the 
possible  results  obtainable  through  the  use  of  this  system. 

The  first  example  is  a  square  (shape  1)  compared  to  a  deformed 
square  (shape  2).  Using  the  line  drawing  capability  to  connect 
points  1  to  2,  2  to  3,  ...,  7  to  8,  and  8  to  1,  the  initial 
shapes  can  be  drawn  using  the  commands  dl_  and  c[2  .  This  is 
shown  in  Figure  1  in  which  numbers  have  been  added  in  order  to 
identify  homologous  points. 

Typing  the  command  ls^  causes  shape  2  to  be  transformed 
to  the  least  squares  fit.  A  display  superimposing  both  shapes 
can  be  obtained  with  the  command  ds  ,  resulting  in  Figure  2a. 
Residual  vectors  starting  from  the  points  of  shape  1  and  ending 
at  the  homologous  points  of  shape  2  can  be  drawn  using  the 
command  d£  .  Figure  2b  shows  these  residuals  superimposed  on 
shape  1,  using  commands  d_r  and  d_l  . 

The  robust  fit  by  repeated  medians  is  produced  by  the 
command  r£  .  Using  commands  ds_  ,  djr  ,  and  d_l_  as  before, 
we  can  obtain  the  displays  of  the  fitted  shapes  and  residuals 
shown  in  Figure  3  for  the  robust  fit. 

Figures  2  and  3  show  that  the  least  squares  and  robust 
fits  can  be  very  different.  The  least  squares  comparison 
indicates  a  complicated  relationship  between  the  two  forms,  with 
differences  of  various  sizes  and  in  various  directions  at  all 
points.  In  contrast,  the  robust  fit  shows  clearly  the  localized 
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FIGURE  2.  Least  squares  fit. 

Shapes  superimposed  (a);  residuals  superimposed  on  shape 
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systematic  nature  of  the  shape  differences.  All  but  two  of  the 
homologous  pairs  are  matched  closely,  and  the  simple  relation¬ 
ship  between  the  shapes  is  apparent. 

Figure  4  shows  the  comparison  of  a  human  skull  to  that  of 
a  chimpanzee.  These  data  are  from  Sneath  (1967)  and  the  example 
is  discussed  in  detail  by  Siegel  and  Benson  (1979).  It  is 
included  here  as  a  demonstration  of  the  graphics  enhancement 
capability.  This  capability  allows  for  outlines  and  other  details 
that  will  not  affect  the  fitting  process  itself  (only  the  homo¬ 
logous  points  are  used  for  that)  but  that  will  be  transformed 
appropriately  in  order  to  produce  a  useful  display. 

No  dashed  lines  are  seen  in  Figure  4  because  the  line  capa¬ 
bility  was  not  used.  However,  both  the  line  and  enhancement 
capabilities  can  be  used  simultaneously  if  desired. 
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4.  IMPLEMENTING  AND  USING  THE  SYSTEM 


To  facilitate  the  use  of  this  program  on  different  computer 
systems,  the  program  is  written  in  standard  FORTRAN  but  four 
subroutines  must  be  supplied  by  the  user:  INPUT,  MOVE,  DRAW, 
and  FINISH.  These  subroutines  are  used  to  link  to  the  user's 
database  and  graphics  capabilities,  which  vary  greatly  from  one 
installation  to  another.  Detailed  requirements  for  each  sub¬ 
routine  ire  included  in  the  program  listing  the  Appendix  after 
each  of  these  four  subroutine  statements. 

The  system  responds  i nteracti vely  to  commands  issued  by  the 
user.  Each  command  is  two  letters  (additional  letters  on  the 
same  line  are  ignored)  and  commands  fall  into  four  groups:  fit, 
draw,  print,  and  set.  Commands  can  be  issued  in  any  order 
desired.  The  fit  and  set  commands  can  be  used  to  transform 
shape  2;  draw  and  print  commands  will  be  based  on  the  most 
recent  orientation  specified  (e.g.  least  squares,  robust,  or 
other).  For  example,  the  sequence  of  commands  £s  dr  d2  rs 
pc  (one  per  line)  will  draw  the  residuals  from  the  least  squares 
fit,  draw  shape  2  (still  in  least  squares  configuration),  and 
finally  print  the  coordinates  of  the  robust  fit. 


5.  COMMAND  DESCRIPTIONS 


There  are  three  fit  commands:  £_s(least  squares),  rs  ( robust ) , 
and  d£(original  data).  Commands  ls_  and  rs  respectively 
produce  the  least  squares  and  robust  fits  discussed  in  Section  2. 
The  original  data  configuration  can  always  be  brought  back  by 
using  da^  .  No  printed  or  graphic  output  results  from  these 

commands,  but  the  internal  configuration  of  shape  2  is  transformed 

accordi ngly . 

There  are  four  draw  commands:  d^(draw  shapes),  dl (draw 
shape  1),  d£(draw  shape  2),  and  dr(draw  residual  vectors). 

Both  shapes  are  drawn  when  ds.  is  used;  shape  1  alone  results 
when  dl_  is  used;  shape  2  alone  is  drawn  when  d2^  is  used. 

There  are  three  components  to  a  shape  drawing:  the  homologous 

points  (crosses  for  shape  1  and  squares  for  shape  2),  straight 
lines  connecting  pairs  of  points  (dashed  lines  for  shape  1, 
solid  lines  for  shape  2),  and  enhancements  (points  and  lines 
that  do  not  affect  the  fitting  process  but  allow  the  inclusion 
of  outlines  and  details  to  be  transformed  appropriately).  The 
residual  vectors  are  arrows  pointing  from  the  points  of  shape  1 
to  the  homologous  points  of  shape  2.  These  vectors  are  drawn 
when  command  d_r  is  given. 

Four  print  commands  are  included:  ££(print  coordinates), 
££(pri nt  parameters ),  £_r  (pri  nt  residuals),  and  h e 1 p ( p r i n t 


command  summary).  The  command  ££  prints  the  coordinates  of 
the  points  of  shape  1,  the  current  values  of  the  transformed 
homologous  points  of  shape  2,  and  the  current  parameter  values. 
The  current  parameter  values  a,b,c,d,p,  and  6,  as  defined 
in  Section  2,  are  printed  when  ££  is  used.  The  command  ££ 
prints  the  residual  vectors  (u(i ) -x ( i )  ,v(i )-y(i ) )  together 
with  their  lengths,  their  average  length,  their  root  mean 
square  length,  and  the  current  parameter  values.  A  summary  of 
valid  commands  is  printed  when  he! p  is  used. 

There  are  four  set  commands:  sjj(set  parameters),  sm 
(modify  parameters),  s£(set  scale),  and  ^e(set  point  and  dash 
size).  These  commands  prompt  you  for  values  to  be  typed  in 
one  per  line  with  five  or  fewer  characters  per  value,  including 
the  required  decimal  point.  You  may  choose  the  transformation 
of  shape  2  directly  by  specifying  the  parameter  values  a,b,p,  an 
using  command  s£  .  To  modify  the  current  parameter  values 
a,b,p,  and  0,  the  command  sm  will  prompt  you  for  changes 
a',b',p',  and  e'  .  The  new  parameter  values  will  then  be 
a+a',  b+b',  pp',  and  0+9'  ,  unless  p's0  in  which  case  p  will 
be  left  unchanged.  The  overall  graphics  size  can  be  specified 
using  the  s_s  command;  the  default  size  of  both  shapes  will 
be  multiplied  by  the  number  entered.  The  size  of  the  homologous 
points  of  both  shapes  and  the  spacing  between  dashes  of  the  lines 
of  shape  1  will  both  be  changed  by  a  factor  entered  using  the 
se  command. 

The  command  en  stops  execution  of  the  program. 
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APPENDI X :  FORTRAN  PROGRAM 


c  interactive  system  for  fitting  and  plotting  shapes, 
c 

c  programmed  lfl/80  by  Andrew  F.  Siegel, 

c 

c  commands  are  listed  by  entering  "help". 

c 

c 

c  variables: 

c  n  =  number  of  h-points  (homologous  points)  of  ®ach  shape, 

c  m  *  number  of  pairs  of  h-points  to  be  connected, 

c  set  m=0  if  no  such  lines  are  desired. 

c  nn  =  500  =  dimension  of  z,  w,  and  ww  (used  for  enhancement 
c  of  drawings) . 

c  ifit  =  1,  2,  3,  or  4  according  to  the  current  fit: 
c  original  data,  least  squares,  robust,  or  special  choice 

c  corns (16)  =  table  of  commands. 

c  sf  ■  scale  factor.  Initially  set  to  one,  can  be  reset  to 
c  change  size  of  drawings, 

c  e  »  point  size  and  dash  spacing. 

c  a,b  *  translation  parameters  to  fit  shape  2  to  shape  1. 
c  rho, theta  *  magnification  and  rotation  parameters 
c  to  fit  shape  2  to  shape  1. 

c  c,d  *  parameters  of  linear  form  of  magnification  and  rotation 
c  to  fit  shape  2  to  shape  1.  (c*rho*cos ( th®ta)  , 

c  d»rho*sin ( theta) )  . 

c  x(n),y(n)  *  coordinates  of  h-points  of  shape  1. 
c  uu(n),vv(n)  *  coordinates  of  original  h-points  of  shape  2. 
c  u(n),v(n)  *  coordinates  of  transformed  h-points  of  shape  2. 
c  lines(m,2)  ■  pairs  of  h-points  to  be  connected  in  drawing: 
c  for  example,  if  lines(k,l)=i  and  lines(k,2)*j 

c  then  the  kth  line  drawn  will  connect  th® 

c  ith  to  the  jth  h-points  of  each  shape, 

c  using  a  solid  line  for  shape  1  and  a  dotted 

c  line  for  shape  2. 

c  z(nn)  *  piecewise  linear  enhancements  for  shape  1.  for 
c  example,  to  draw  an  outline  connecting  (al,bl)  to 

c  (a2,b2)  to  ...  to  (ak,bk),  store  k .al ,bl , . . . ,ak,bk 

c  in  the  vector  z() .  Additional  outlines  can  be 

c  added  one  after  another  in  this  same  form,  endinq 

c  finally  with  a  zero. 

c  ww(nn)  »  original  piecewise  linear  enhancements  for  shape  2. 

c  w(nn)  *  transformed  piecewise  linear  enhancements  for  shape  2 
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c  ut(n),vt(n)  =  temporary  storage  for  repeated  median  fitting. 

c 

c 

dimension  x(100) ,y(100) ,u{100) ,v(100) ,uu(100) ,vv(100) ,1 ines ( 200 , 2) 

*  ,2(500) , w ( 500 ) , ww ( 5  0  0 ) ,coms (16 ) ,ut(100) ,vt(100) 
data  corns/ 'Is' ,  ’  rs ' ,'da' ,'ds' ,'dl' ,'d2' ,'dr' ,’pc' ,'pp' ,'pr' , '  he ' 

*  , '  se '  , 1  sp '  , '  sm '  , '  ss  '  , '  en '  / 
write(6,l) 

1  format ('  enter  11  help'1  for  a  list  of  commands.’) 
nn*500 
z(l)*0. 
ww( 1 ) *0 . 
sf-1. 
e* . 0075 

call  input (n ,nn ,m,x ,y ,uu ,vv , lines ,z ,ww) 
if it»l 
a*0 . 
b»0. 
c*l . 
d«0 . 

call  coord (n,nn,a,b,c,d,u,v ,uu , vv ,w, ww) 

10  read(5,2)com 

2  format(a2) 

if (com .eq. corns (1))  call  ls(n,nn,a,b,cfd,ifit,x,y,u,v,uu,vv ,w,ww) 
if (com. eq. corns (2) )  call  rs (n ,nn ,a ,b ,c ,d , if it ,x ,y ,u ,v,uu,vv ,w,ww, 

*  ut,vt) 
if (com.eq.coms(3) )  call  data (n ,nn ,a ,b,c ,d , if it , u ,v ,uu ,vv ,w,ww) 
if (com .eq. corns (4) )  call  ds (0 ,n ,nn ,m, if it ,sf ,e ,x ,y,u ,v, lines ,w, z) 
if ( com .eq. corns (5))  call  ds(l,n,nn,m,ifit,sf,e,x ,y ,u ,v , 1 ines ,w,z) 
if ( com. eq. corns (6))  call  ds(2,n,nn,m,ifit,sf,e,x ,y ,u ,v ,1 ines ,w, z) 
if (com.eq.coms(7) )  call  dr (n ,nn , if  it , sf ,e ,x ,y ,u  ,v  ,z) 
if (com. eq. corns (8) )  call  pc (n ,a ,b ,c ,d , if  it ,x ,y ,u ,v) 
if (com.eq.coms(9) )  call  pp(a,b,c,d) 
if (com.eg.coms(10) )  call  pr (n ,a ,b , c ,d , if it ,x ,y,u,v) 
if (com.eq.coms(ll) )  call  help 
if (com. eq. corns (12) )  call  se(e) 

if (com.eq.coms(13) )  call  sp (n ,nn ,a ,b ,c ,d , if it ,u ,v ,uu ,vv , 

*  w,ww) 
if  (com  .eq. corns  (14) )  call  sm  (n  ,nn  ,a  ,b  ,c  ,d ,  if  it  ,u  ,v,uu  ,vv , 

*  w,ww) 


no  ooooooo  no  oooooono  o  on 
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if (com.eq.coms(15) )  call  ss(sf) 

if (com. eq. corns (16) )  stop 

go  to  10 

stop 

end 


subroutine  input (n ,nn ,m ,x ,y ,uu ,vv , 1 ines ,z ,ww) 
reads  in  coordinate  and  enhancement  data. 

dimension  x(n) ,y(n) ,uu(n) ,vv(n) ,lines(m,2) ,z(nn) ,ww(nn) 

***************************************************** 
the  user  must  supply  this  subroutine  to  obtain  the 
number  of  points,  n;  the  initial  coordinates  x(),y(), 
uu()  ,  and  vv();  the  values  of  m  and  lines(,);  and 
the  enhancement  values  z()  and  ww(). 
***************************************************** 


return 

end 


subroutine  move(x,y) 
moves  cursor  to  (x,y) . 

***************************************************** 

the  user  must  supply  this  subroutine  to  move  the 
graphics  cursor  to  (x,y)  but  draw  nothing,  generally, 
x  and  y  will  lie  between  0  and  1,  although  these 
limits  can  be  exceeded. 

***************************************************** 

return 

end 


subroutine  draw(x,y) 
c  draws  a  line  to  (x,y) . 

c  ***************************************************** 
c  the  user  must  supply  this  graphics  subroutine  to  draw 

c  a  line  from  the  previous  cursor  location  to  the  point 

c  (x,y) .  the  cursor  location  should  then  be  updated 


n o  ooooo  o  o  on 
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to  (x,y) . 

***************************************************** 

return 

end 


subroutine  finish 
turns  off  graphics  mode. 

***************************************************** 

the  user  must  supply  this  subroutine  to  turn  off 
graphics  mode,  if  needed. 

***************************************************** 

return 

end 


subroutine  Is  (n  ,nn  ,a  ,b,c  ,d  ,  if  it  ,x  ,y ,  u  ,v  ,uu  ,vv  ,w,  ww) 
dimension  x{n) ,y{n)  ,u(n)  ,v(n)  ,uu(n)  ,vv(n)  ,w(nn) ,ww(nn) 
if it*2 
an»n 
sx*0 . 
sy-0 . 
su*0. 
sv«0 . 
sux»0 . 
suy-0 . 
suu*0. 
svx*0. 
svy«0 . 
svv*0 . 
do  10  i-l,n 
sx»sx+x ( i) 
sy*sy+y ( i) 
su»su+uu { i) 
sv»sv+vv { i) 
sux»sux+uu(  i)  *x(  i) 
suy«suy+uu ( i) *y( i) 
suu«suu+uu { i) *uu ( i) 
svx»svx+w  ( i)  *x(  i) 
svy»svy+vv ( i)  *y  ( i) 
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10  svv*svv+vv { i) *vv ( i) 

ss*suu-su**2/an+svv-sv**2/an 
c* (sux-su*sx/an+svy-sv*sy/an) /ss 
d* (suy-su*sy/an-svx+sv*sx/an) /ss 
a* (sx-c*su+d*sv) /an 
b* (sy-c*sv-d*su) /an 

call  coord (n  ,nn  ,a  ,b , c  ,d  ,u  ,v  ,uu  ,vv  ,w,ww) 

return 

end 


subroutine  rs(n,nn,a,b,c,d,ifit,x,y,u,v,uu,vv,w,ww,ut,vt) 
c  robust  fit  by  repeated  medians. 

dimension  x(n) ,y(n) ,u(n) ,v(n) ,uu(n) ,vv(n) ,w(nn) ,ww(nn) ,ut(n) ,vt(n) 
call  lstnjnnja/bfCjdjifitjX^jUfVjUUfWjW/Ww) 
if  it*  3 

c  get  magnification  rho 
do  20  i*l,n 
jl*0 

do  10  j*l,n 

if(i.eq.j)  go  to  10 
jl*jl+l 

ut( jl)*sqrt( ( (x{ j)-x(i) ) **2+(y( j) -y(i) ) **2 ) / 

*  ((u(j)-u(i))**2+(v(j)-v(i))**2)) 

10  continue 

20  vt( i) *fmed (n-1 ,ut) 

rho*fmed  (n,vt) 
c  get  rotation  theta 
do  40  i*l,n 
jl*0 

do  30  j*l,n 

if(i.eq.j)  go  to  30 
jl-jl+1 
xl*x ( j ) -x ( i) 
yl*y{ j)-y(i) 
ul*u( j)-u(i) 
vl*v( j)-v(i) 

ut( jl)-atan2(ul*yl-vl*xl,ul*xl+vl*yl) 
continue 

vt( i) *fmed (n-1 ,ut) 


30 

40 
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theta*fmed (n ,vt) 
c  get  translation  vector  a,b 
ct»rho*cos (theta) 
dt»rho*sin ( theta) 
do  50  i*l,n 

ut  (  i)  »x  (  i)  -ct*u(  i)  +dt*v  (  i) 

50  vt(  i)  »y(  i)  -dt*u(  i)  -ct*v(  i) 

at*fmed (n  ,ut) 
bt*fmed  (n  ,vt) 
at*at+ct*a-dt*b 
b»bt+dt*a+ct*b 
a=at 

ctt»ct*c-dt*d 

d«ct*d+dt*c 

c«ctt 

call  coord(n,nn,a,b,c,d,u,v,uu,vv,w,ww) 

return 

end 


subroutine  data ( n , nn ,a ,b, c ,d,ifit,u,v,uu ,vv,w,ww) 
puts  original  data  back  for  shape  2. 

dimension  uu{n) ,vv(n) ,u(n) ,v(n) ,ww(nn) ,w(nn) 

ifit-1 

a«0 . 

b»0 . 

c»l . 

d*0 . 

call  coord (n ,nn ,a ,b,c ,d ,u ,v ,uu ,vv ,w,ww) 

return 

end 


subroutine  coord(n,nn,a ,b,c,d,u,v,uu,vv,w,ww) 
c  puts  transformed  coordinates  into  u,v,w  for  shape  2. 
dimension  u(n) ,v(n) ,uu(n) ,vv(n) ,w(nn) ,ww(nn) 
do  10  i«l,n 

u(  i)  «a+c*uu  (  i)  -d*w  ( i) 
v ( i) »b+d*uu ( i) +c*vv ( i) 
i«l 


10 
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20  j*ww(i) 
w( i) »j 
i*  i+1 

if(j.eq.0)  return 
do  30  k*l ,  j 

w(  i)  *a+c*ww(  i) -d*ww( i+1) 
w( i+1) =b+d*ww( i) +c*ww( i+1) 
30  i-i+2 

go  to  20 
end 


subroutine  ds (nshape ,n,nn,m,ifit,sf,e,x ,y ,  u ,v ,1 ines ,w, z) 
c  draw  shapes  (both  if  nshape=0,  shape  1  only  if  nshape=l, 
c  and  shape  2  only  if  nshape*2) . 

dimension  x(n) ,y(n)  fu(n) ,v(n) ,lines(m,2) ,w(nn) ,z(nn) 
call  desc(ifit) 
call  scale(n,s,sf ,xm,ym,x,y) 
if (nshape. eq. 2)  go  to  40 
c  draw  shape  1 

do  10  i*l,n 

p* (x(i) -xm)/s+.5 
q»(y(i)-ym)/s+.5 
call  move(p,q+e) 
call  draw(p,q-e) 
call  move(p+e,q) 

10  call  draw(p-e,q) 

call  enh(nn,xm,ym,sfz) 
if(m.eq.0)  go  to  40 
do  30  i*l,m 
il«lines ( i ,1) 
i2*lines ( i ,2) 
xl*(x(il) -xm)/s+. 5 
yl*(y( il) -ym) /s+. 5 
x2*(x(i2)-xm)/s+.5 
y2*(y(i2)-ym)/s+. 5 
s2«sqrt ( (x2-xl) **2  + (y2-yl) **2) 
if (s2.lt. 4 .*e)  go  to  30 
j* (s2-3 . *e) /( 2 . *e) 
xu« (x2-xl)/s2 
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yu* (y2-yl)/s2 

exu«e*xu 

eyu*e*yu 

exu2*2.*exu 

eyu2»2.*eyu 

p»(s2-3.*e-2.*e*float{ j) )/2. 

xl*xl+p*xu+exu2 
yl»yl+p*yu+eyu2 
do  20  k*l,j 

call  move(xl,yl) 
call  draw(xl+exu,yl+eyu) 
xl*xl+exu2 
20  yl=yl+eyu2 

30  continue 

40  if (nshape.eq.l)  go  to  70 
c  draw  shape  2 

do  50  i»l»n 

p*(u(i)-xm)/s+.5 
q* (v( i) -ym)/s+. 5 
call  move(p+e,q+e) 
call  draw(p+e ,q-e) 
call  draw(p-e ,q-e) 
call  draw(p-e ,q+e) 

50  call  draw{p+e,q+e) 

call  enh(nn,xm,ym,s,w) 
if (m.eq.0)  go  to  70 
do  60  i*l,m 
il»lines( i ,1) 
i2*lines ( i , 2) 
xl«(u(il)-xm)/s+.5 
yl* (v( il) -yro) /s+ . 5 
x2*(u(i2)-xm)/s+.5 
y2*(v(i2)-ym)/s+.5 
s2»sqrt< (x2-xl)**2+(y2-yl)**2) 
if (s2.lt. 4. *e)  go  to  60 
xu* (x2-xl) /s2 
yu* (y2-yl) /s2 

call  move (xl+2 .*e*xu,yl+2 .*e*yu) 
call  draw(x2-2.*e*xu,y2-2.*e*yu) 
60  continue 
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70  call  finish 

return 
end 
c 
c 

subroutine  enh (nn ,xm ,ym  ,s  , z) 
c  draws  enhancements 
dimension  z(nn) 
i*l 

10  j*z(i) 

i=i+l 

if(j.eq.0)  return 
xl*(z(i)-xm)/s+.5 
yl* (z ( i+1) -ym) /s+.  5 
call  move{xl,yl) 
call  draw(xl,yl) 
i=i+2 

if(j.eq.l)  go  to  10 
do  20  k-2,j 

call  draw ( ( z ( i) -xm) /s+. 5 , ( z ( i+1) -ym) /s+. 5) 
20  i*i+2 

go  to  10 
end 


subroutine  dr(n,nn,ifit,sffe,x,y,u,v,z) 
c  draws  residual  vectors 

dimension  x(n)  ,y(n)  ,u(n)  ,v(n)  ,z(nn) 
call  desc(ifit) 
call  scale(n,s,sf ,xm,ym,x,y) 
do  10  i*l,n 

xl*  (x ( i) -xm) /s+ . 5 
yl*(y(i)-ym)/s+.5 
ul*(u( i) -xm)/s+. 5 
vl«(v(i) -ym)/s+.5 
call  move(xl,yl) 
call  draw(ul,vl) 
dx« (ul-xl) /8 . 
dy- (vl-yl) /8 . 

call  draw(xl+7 . *dx+dy ,yl+7, *dy-dx) 
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call  move(ul,vl) 

10  call  draw(xl+7 . *dx-dy,yl+7 . *dy+dx) 

call  finish 
return 
end 


subroutine  pc(n,a,b,c,d,ifit,x,y,u,v) 
prints  coordinates  and  parameters, 
dimension  x(n) ,y(n) ,u(n) ,v(n) 
call  desc(ifit) 

write (6,1)  (x(i),y(i),u(i),v(i) , i*l,n) 

format ('  coordinates  x,y  of  shape  1  and  u,v  of  shape 
*  ( lx , 4f 10 . 5 ) ) 

call  pp(a,b,c,d) 
return 
end 


subroutine  pr(n,a,b,c,d,ifit,x,y,u,v) 

c  print  residuals 

dimension  x(n) ,y(n) ,u(n) ,v(n) 
call  desc(ifit) 
write(6 ,1) 

1  formate  residuals:  coordinates  and  lengths.') 
sr*0 . 

srr*0 . 
do  10  i»l,n 
rx*u( i) -x ( i) 
ry*v  (  i)  -y  (  i) 
r«sqrt ( rx**2+ry**2) 
write(6,2)  rx,ry,r 

2  format(lx,3fl0 .5) 
sr»sr+r 

10  srr»srr+r**2 

sr*sr/f loat (n) 
srr*sqrt(srr/float(n) ) 
write(6 ,3) sr  ,srr 

3  format(/'  average  residual  length  *  ',fl0.5/ 

*  '  r.m.s.  residual  length  ■  ' rfl0.5) 
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call  pp(a,b,c,d) 

return 

end 


subroutine  pp(a,b,c,d) 
print  parameters. 

rho*sqrt (c**2+d**2) 
theta*atan2 (d,c) *180 ./3 . 141593 
wr ite (6 ,1) a  ,b  ,c ,d ,rho  , theta 

formate  parameters  a,  b,  c,  d,  rho,  theta  (in  degrees) 
*  /lx , 6f 10 . 5 ) 

return 
end 


subroutine  help 
write (6,1) 

1  format ('  valid  commands  are:'// 

*  5x, ' fit: */10x, 'least  squares ', 17x ,' Is ' / 

*  10x, ’robust  (repeated  median) ' ,6x, *rs'/ 

*  10x, 'original  data’ ,17x, 'da'// 

*  5x, 'draw: '/10x, 'shapes' ,24x, 'ds'/ 

*  10x, 'shape  1  only' ,16x , ' dl ' / 

*  10x, 'shape  2  only' ,16x, 'd2 '/ 

*  10x, ' residual  vectors ', 14x ,' dr ' // 

*  5x,' print: '/10x, 'coordinates' ,19x,'pc'/ 

*  10x ,' parameters ', 20x , 1 pp' / 

*  10x, ' residuals' ,21x ,' pr '/ 

*  10x, 'commands' ,20x, 'help'// 

*  5x, 'set: '/10x, 'parameters  (new) ' ,14x, ' sp’/ 

*  10x, ' parameters  (mod ify) ' , llx , ’ sm ' / 

*  10x, 'scale' ,25x,’ss'/ 

*  10x, 'point  size,  dash  spacing’ ,6x, 'se'// 

*  5x,  'end' ,32x, 'en'/) 
return 

end 


c 

c 


subroutine  se(e) 
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c  set  point  size-  and  dash  spacing  for  drawing  shapes, 
write (6 , 1) 

1  format^'  please  enter  desired  dash  and  point  size  for  drawing.’) 
read { 5 , 2) e 

2  format ( f 5 . 2) 
e* . 0075*e 
write (6, 3) 

3  format ( '  thank  you.') 
return 

end 

c 

c 

subroutine  ss(sf) 

c  sets  the  scale  factor  to  any  desired  magnification  for  drawing, 
write (6 ,1) 

1  format('  please  enter  desired  magnification  factor' 

*  , 1  for  drawing . ' ) 
read (5 ,  2)  sf 

2  format(f5.2) 
write (6, 3) 

3  format ('  thank  you.’) 
return 

end 

c 

c 

subroutine  sp(n,nn,a,b,c,d,ifit,u,v,uu  rvv ,w,ww) 
c  sets  transformation  parameters  to  any  desired  values, 
dimension  u(n) ,v(n) ,uu(n) ,vv(n) ,w(nn) ,ww(nn) 
if it*4 
write (6 ,1) 

1  format('  desired  transformation  parameters...’/ 

*  '  enter  a,  b,  rho,  and  theta  (degrees),  one  per  line') 
read (5,2)a,b,rho,  theta 

2  format(f5.2) 

if (rho.eq.0. )  rho*l. 
theta* theta *3 . 141593/180 . 
c* rho* cos ( theta) 
d*rho*sin ( theta) 

call  coord  (n  ,nn  ,a  ,b  ,c  ,d  , u  ,  v  ,uu  ,vv  ,w  ,ww) 
write(6,3) 
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3  formate  thank  you.') 

return 
end 
c 
c 

subroutine  sm(n,nn,a,b,c,d,ifit,u,v,uu,vv , w ,ww) 
c  modify  current  values  of  transformation  parameters, 
dimension  u(n)  ,v{n) ,uu(n)  fvv(n)  ,w(nn)  ,ww(nn) 
if it=4 
write(6,l) 

1  format('  modification  of  current  transformation  parameters...' 

*  /'  enter  a,  b,  rho,  and  theta  (degrees)  ,  one  per  line, 

read (5 , 2) at ,bt , rho , theta 

2  format ( f 5 . 2) 

if (rho .eq.0. )  rho=l. 

theta* theta *3 .141593 /a 80. 

ct*rho*cos ( theta) 

dt*rho*sin ( theta) 

att=at+a*ct-b*dt 

b*bt+a*dt+b*ct 

a*att 

ctt=c*ct-d*dt 

d*c*dt+d*ct 

c=ctt 

call  coord (n ,nn ,a ,b ,c ,d , u ,v ,uu ,vv ,w,ww) 
write(6 ,3) 

3  formatC  thank  you.') 
return 

end 


subroutine  desc(ifit) 
c  prints  the  current  type  of  fit. 

go  to  (10,20,30,40) ifit 
10  write(6,l) 

1  format ('  original  data') 

return 

20  write(6,2) 

2  formate  least  squares  fit') 

return 
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30  write (6 , 3 ) 

3  format('  robust  fit  by  repeated  medians') 

return 

40  write(6,4) 

4  formate  special  fit') 
return 

end 


subroutine  scale (n  ,s , sf ,xm ,ym ,x ,y) 
finds  center  and  scale  for  drawing  so  that  shape  1  fits  in  square 
from  (.25,. 25)  to  (.75, .75),  then  adjusts  by  the  value  of  sf. 
dimension  x(n),y(n) 
xmin*x(l) 
xmax=x (1) 
ymin=y(l) 
ymax=y(l) 
do  10  i*l,n 

xmin*aminl  (xmin  ,x  ( i)  ) 
xmax=amaxl (xmax ,x ( i) ) 
ymin*aminl  (ymin  ,y  (  i) ) 

0  ymax^amaxl (ymax ,y ( i) ) 

s*2. *amaxl (xmax-xmin ,ymax-ymin) /sf 

xm= (xmin+xmax) /2 . 

ym=*  (ymin+ymax)  /  2 . 

return 

end 


function  fmed(n,a) 

c  finds  the  median  of  a(l),...,a(n). 
c  a  more  efficient  algorithm  could  be  used 
c  if  n  is  large  to  decrease  execution  time, 
dimension  a(n) 
m«n/2+l 
do  20  i»l,m 
k«i 

il»i+l 

do  10  j*il,n 

if (a( j) .lt.a(k) )  k»j 
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at*a  { i) 
a(  i)  *a  (  k) 

20  a(k)-at 

if (n.eq.2*(n/2) )  go  to  30 

fmed*a (m) 

return 

30  fmed« (a(m) +a (m-1) ) /2 . 

return 
end 
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